###################################################################################################
## Payson (2021) When Cities Lobby Replication Code
## Chapter Three: Figures (3.1 - 3.8) 
## Created with R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
###################################################################################################


## Load packages
library(ggplot2); library(RColorBrewer); library(effects); 
library(tidyverse); library(ggpubr)


## Functions for clean tex tables (need to define for tables to compile)
pretty <- function(x) {
  paste(prettyNum(x, big.mark = ","))}

num.cities <- function(x) {
  length(which(getfe(x)=="fips"))}

N <- function(x) {
  x$N}

getmean <- function(x) {
  round(mean(x$response), 2)}


## Set working directory
setwd()


###################################################################################################
## Figure 3.1 Correlates  of  Lobbying  Across  Cities
###################################################################################################

dta <- read.csv("city_lobby.csv")


## Specify model to get marginal effects
m <- lm(lobby ~ lower.num.reps + np.distance + log.pop + pct.white + 
          log.income + I(log.income^2) + log.own.source.pp + 
          log.med.house + factor(year) + state, data = dta)


## Get predicted values from model
pred <- predict(m)


## Merged predicted values with covariate values
plot.dta <- data.frame(pred = pred, log.pop = m$model$log.pop, 
                       num.reps = m$model$lower.num.reps,
                       np.distance = m$model$np.distance,
                       med.house = m$model$log.med.house,
                       own.source = m$model$log.own.source.pp,
                       income = m$model$log.income,
                       pct.white = m$model$pct.white)



## Predicted probability of lobbying by population
(p1 <- ggplot(plot.dta, aes(x = log.pop, y = pred)) +
    #  geom_point(alpha = .1) +
    stat_summary_bin(geom = "point", bins = 200) +
    stat_smooth(method = "lm", formula = y ~ poly(x, 3), 
                color = "black", se = FALSE) +
    ylim(0, 1) + xlim(9, 15) +
    labs(y = "Probability of Lobbying", x = "Population (Log)")  +
    theme_bw(base_family = "serif", base_size = 11) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))


## Predicted probability of lobbying by own source revenue per capita
(p2 <- ggplot(plot.dta, aes(x = own.source, y = pred)) +
    #  geom_point(alpha = .01) +
    stat_summary_bin(geom = "point", bins = 100) +
    stat_smooth(method = "lm", formula = y ~ poly(x, 2), 
                color = "black", se = FALSE) +
    coord_cartesian(ylim = c(0, .3)) + xlim(6, 9) +
    labs(y = "Probability of Lobbying", x = "Own Source Revenue Per Capita (Log)")  +
    theme_bw(base_family = "serif", base_size = 11) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))


## Predicted probability of lobbying by income
(p3 <- ggplot(plot.dta, aes(x = income, y = pred)) +
    stat_summary_bin(geom = "point", bins = 100) +
    geom_smooth(method = lm, formula = y ~ poly(x, 2), 
                color = "black", se = FALSE) +
    coord_cartesian(ylim=c(0.0,.4)) + xlim(10, 12) +
    labs(y = "Probability of Lobbying", x = "Income (Log)")  +
    theme_bw(base_family = "serif", base_size = 11) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))



## Predicted probability of lobbying by median home value
(p4 <- ggplot(plot.dta, aes(x = med.house, y = pred)) +
    stat_summary_bin(geom = "point", bins = 100) +
    stat_smooth(method = "lm", formula = y ~ poly(x, 2), 
                color = "black", se = FALSE) +
    xlim(11, 13.5) +
    labs(y = "Probability of Lobbying", x = "Median House Value (Log)")  +
    theme_bw(base_family = "serif", base_size = 11) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))


## Combine plots
(fig <- ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2))



###################################################################################################
## Figure 3.2 Representational  Correlates  of  Lobbying
###################################################################################################

dta <- read.csv("city_lobby.csv")


## Specify model to get marginal effects
m <- lm(lobby ~ lower.num.reps + np.distance + log.pop + pct.white + 
          log.income + I(log.income^2) + log.own.source.pp + 
          log.med.house + factor(year) + state, data = dta)


## Get predicted values from model
pred <- predict(m)


## Predicted probability of lobbying by number of representatives
(p1 <- ggplot(plot.dta, aes(x = num.reps, y = pred)) +
   xlim(0, 30) + ylim(0, 1) +
   stat_summary_bin(geom = "point", bins = 30) +
   geom_smooth(method = lm, formula = y ~ poly(x, 2), 
               color = "black", se = FALSE) +
   labs(y = "Probability of Lobbying", x = "Number of Lower House Reps")  +
   theme_bw(base_family = "serif", base_size = 11) +
   theme(plot.background = element_blank(),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank()))


## Predicted probability of lobbying by ideological distance between reps
(p2 <- ggplot(plot.dta, aes(x = np.distance, y = pred)) +
    stat_summary_bin(geom = "point", bins = 10) +
    geom_smooth(method = lm, formula = y ~ poly(x, 2), 
                color = "black",se = FALSE) + xlim(0, 4) +
    labs(y = "Probability of Lobbying", x = "Ideological Distance Between Reps")  +
    theme_bw(base_family = "serif", base_size = 11) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))

(fig <- ggarrange(p1, p2, ncol = 2, nrow = 1))



###################################################################################################
## Figure 3.3 Effect of Gaining More Representatives After Redistricting on Lobbying
###################################################################################################

dta <- read.csv("city_lobby.csv")

## For each year, get mean lobbying probability for cities that gained reps through redistricting
## vs. cities that didn't gain reps through redistricting 

tmp <- dta %>%
  group_by(year, redistricted) %>%
  summarize(lobby = mean(lobby, na.rm = T))


## Subset to complete data
tmp <- tmp[tmp$year>2005,]

tmp$redistricted <- ifelse(tmp$redistricted==1, "More Reps After Redistricting", "Same Number of Reps")

(p <- ggplot(tmp, aes(year, lobby, linetype = factor(redistricted))) +
    geom_point() +
    geom_line() +
    geom_vline(xintercept = 2010) +
    scale_x_continuous(breaks = c(2006, 2008, 2010, 2012, 2014)) +
    scale_y_continuous(breaks = c(.14, .16, .18, .20, .22, .24, .26)) +
    theme_bw(base_family = "serif", base_size = 12) +
    theme(panel.spacing = unit(2, "lines"),
          legend.title = element_blank(),
          legend.position = "top",
          plot.margin = unit(c(1,2,1,1), "cm"),
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white")) +
    labs(y = "Probability of Lobbying", x = ""))



###################################################################################################
## Figure 3.4 Classifying Republican and Democratic Cities
###################################################################################################

dta <- read.csv("mrp.csv")


year2008 <- subset(dta, year == 2008)

NY <- year2008[year2008$city == "NEW YORK CITY", ]
Berkeley <- year2008[year2008$city == "BERKELEY CITY" & year2008$state=="California", ]
Amarillo <- year2008[year2008$city == "AMARILLO CITY", ]
Mesa <- year2008[year2008$city == "MESA CITY", ]
Jacksonville <- year2008[year2008$city == "JACKSONVILLE CITY" & year2008$state=="Florida", ]
Auburn <- year2008[year2008$city == "AUBURN CITY" & year2008$state=="Alabama", ]
Iowa <- year2008[year2008$city == "IOWA CITY CITY", ]


(p <- ggplot(year2008, aes(x = pres.dem.2008, y = (-mrp.estimate), color = both)) +
    geom_point(alpha = .7) +
    stat_smooth(method = "lm", color = "black") +
    labs(y = "City Liberalism Score", x = "Democratic Vote Share (2008)") +
    theme_bw(base_family = "serif", base_size = 12) +
    scale_colour_grey(name="", start=0.02, end = .8, na.translate = FALSE) +
    theme(legend.position = "top",
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()) +
    guides(color = guide_legend(reverse = TRUE)) +
    annotate("text", x = NY$pres.dem.2008, y = -(NY$mrp.estimate) + .1, 
             label = "New York City", size = 4) +
    annotate("text", x = Berkeley$pres.dem.2008 - .03, y = -(Berkeley$mrp.estimate) + .1, 
             label = "Berkeley, CA", size = 4) +
    annotate("text", x = Amarillo$pres.dem.2008 + .05, y = -(Amarillo$mrp.estimate) - .05, 
             label = "Amarillo, TX", size = 4) +
    annotate("text", x = Mesa$pres.dem.2008, y = -(Mesa$mrp.estimate), 
             label = "Mesa, AZ", size = 4) +
    annotate("text", x = Jacksonville$pres.dem.2008 + .02, y = -(Jacksonville$mrp.estimate), 
             label = "Jacksonville, FL", size = 4) + 
    annotate("text", x = Auburn$pres.dem.2008, y = -(Auburn$mrp.estimate), 
             label = "Auburn, AL", size = 4) +
    annotate("text", x = Iowa$pres.dem.2008 -.02, y = -(Iowa$mrp.estimate) + .05, 
             label = "Iowa City", size = 4))



###################################################################################################
## Figure 3.5: Effect of Delegation Mismatch on City Lobbying: Difference-in-Differences Predictions
###################################################################################################

dta <- read.csv("mrp.csv")


## Create factors for year and delegation mismatch treatment
dta$year <- factor(dta$year)
dta$mismatch <- factor(dta$mismatch)


## Estimate models for full sample and by city size
m <- lm(lobby ~ mismatch + log.pop + log.income + log.own.source.pp + log.med.house
        + fips + year, data = dta)

m1 <- lm(lobby ~ mismatch + log.pop + log.income + log.own.source.pp + log.med.house
         + fips + year, data = subset(dta, pop.cat == "medium"))

m2 <- lm(lobby ~ mismatch + log.pop + log.income + log.own.source.pp + log.med.house
         + fips + year, data = subset(dta, pop.cat == "big"))

m3 <- lm(lobby ~ mismatch + log.pop + log.income + log.own.source.pp + log.med.house
         + fips + year, data = subset(dta, pop.cat == "urban"))


## Use "effects" package to estimate marginal effect of lobbying by city size
plot.dta <- data.frame(citysize = c(rep("Full Sample", 2), rep("20,000 - 50,000", 2), 
                                    rep("50,000 - 100,000", 2), rep("Over 100,000", 2)),
                       lower = c(summary(effect("mismatch", m))$lower,
                                 summary(effect("mismatch", m1))$lower,
                                 summary(effect("mismatch", m2))$lower,
                                 summary(effect("mismatch", m3))$lower),
                       upper = c(summary(effect("mismatch", m))$upper,
                                 summary(effect("mismatch", m1))$upper,
                                 summary(effect("mismatch", m2))$upper,
                                 summary(effect("mismatch", m3))$upper),
                       effect = c(summary(effect("mismatch", m))$effect,
                                  summary(effect("mismatch", m1))$effect,
                                  summary(effect("mismatch", m2))$effect,
                                  summary(effect("mismatch", m3))$effect),
                       mismatch = rep(c("Aligned", "Mismatched"), 4))


## Set up data for plotting
plot.dta$citysize <- factor(plot.dta$citysize, levels = c("Full Sample", "20,000 - 50,000", 
                                                          "50,000 - 100,000", "Over 100,000"))
plot.dta$mismatch <- factor(plot.dta$mismatch, levels = c("Aligned", "Mismatched"))
dodge <- position_dodge(width = 0.5)


(p <- ggplot(plot.dta, aes(citysize, effect, fill = mismatch)) +
    geom_bar(position = "dodge", stat = "identity", width = .5) +
    geom_errorbar(aes(ymin = lower, ymax = upper), position = dodge, width = 0.1) +
    scale_fill_grey(name = "", start = .8, end = .2) +
    theme_bw(base_family = "serif",
             base_size = 12) +
    theme(legend.position = "top",
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white")) +
    labs(y = "Probability of Lobbying", x = "City Size"))



###################################################################################################
## Figure 3.6 Effect of Delegation Mismatch on Lobbying by City Partisanship: 
## Difference-in-Differences Predictions
###################################################################################################

dta <- read.csv("mrp.csv")

## Create factors for year and delegation mismatch treatment
dta$year <- factor(dta$year)
dta$mismatch <- factor(dta$mismatch)

## Estimate models by city partisanship
m <- lm(lobby ~ mismatch + log.pop + log.income + log.own.source.pp + log.med.house
        + fips + year, data = subset(dta, both == "Republican City"))

m1 <- lm(lobby ~ mismatch + log.pop + log.income + log.own.source.pp + log.med.house
         + fips + year, data = subset(dta, both == "Democratic City"))


## Use "effects" package to estimate marginal effect of lobbying by city partisanship
plot.dta <- data.frame(type = c(rep("Republican City", 2), rep("Democratic City", 2)),
                       lower = c(summary(effect("mismatch", m))$lower,
                                 summary(effect("mismatch", m1))$lower),
                       upper = c(summary(effect("mismatch", m))$upper,
                                 summary(effect("mismatch", m1))$upper),
                       effect = c(summary(effect("mismatch", m))$effect,
                                  summary(effect("mismatch", m1))$effect),
                       mismatch = rep(c("Aligned", "Mismatched"), 2))


dodge <- position_dodge(width = 0.5)

## Figure 4: Lobbying and State Transfers by City Revenue Tercile
(p <- ggplot(plot.dta, aes(type, effect, fill = mismatch)) +
    geom_bar(position = "dodge", stat = "identity", width = .5) +
    geom_errorbar(aes(ymin = lower, ymax = upper), position = dodge, width = 0.1) +
    scale_fill_grey(name = "", start = .8, end = .2) +
    theme_bw(base_family = "serif", base_size = 12) +
    theme(legend.position = "top",
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white")) +
    labs(y = "Probability of Lobbying", x = "City Partisanship"))



###################################################################################################
## Figure 3.7 Delegation Ideology by Aligned and Mismatched Cities
###################################################################################################

dta <- read.csv("mrp.csv")

## Create factor for delegation mismatch treatment
dta$mismatch <- factor(dta$mismatch)

## Subset to Republican and Democratic cities with no missing values
plot.dta <- subset(dta, both != "Moderate City" & !is.na(lower.mean.np) &
                     !is.na(mismatch))

(p <- ggplot(plot.dta, aes(lower.mean.np, linetype = mismatch)) +
    geom_density(alpha = .8) +
    facet_wrap(~both) +
    scale_linetype_manual(values=c("solid", "dashed"), name = "", na.translate = FALSE, 
                          labels = c("Aligned", "Mismatched")) +
    labs(y = "Density", x = "Average Delegation Conservatism Score") +
    theme_bw(base_family = "serif", base_size = 12) +
    theme(legend.position = "top",
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white")))



###################################################################################################
## Figure 3.8 Figure 3.8 Effect of Delegation Conservatism on Lobbying: Difference-in-Differences 
## Predictions
###################################################################################################

dta <- read.csv("mrp.csv")


## Create factors for year and delegation mismatch treatment
dta$year <- factor(dta$year)
dta$mismatch <- factor(dta$mismatch)


## Estimate models
m <- lm(lobby ~ lower.mean.np + mismatch + log.pop + log.income + log.own.source.pp + log.med.house
        + fips + year, data = subset(dta, both == "Republican City"))


m1 <- lm(lobby ~ lower.mean.np + mismatch + log.pop + log.income + log.own.source.pp + log.med.house
         + fips + year, data = subset(dta, both == "Democratic City"))

np <- list(lower.mean.np=c(-2, -1, 0, 1, 2))


## Use "effects" package to estimate marginal effect of lobbying by city partisanship
plot.dta <- data.frame(rep.type = rep(c(-2, -1, 0, 1, 2), 2),
                       city.type = c(rep("Republican City", 5), rep("Democratic City", 5)),
                       lower = c(summary(effect("lower.mean.np", m, xlevels = np))$lower,
                                 summary(effect("lower.mean.np", m1, xlevels = np))$lower),
                       upper = c(summary(effect("lower.mean.np", m, xlevels = np))$upper,
                                 summary(effect("lower.mean.np", m1, xlevels = np))$upper),
                       effect = c(summary(effect("lower.mean.np", m, xlevels = np))$effect,
                                  summary(effect("lower.mean.np", m1, xlevels = np))$effect))


limits <- aes(ymin = lower, ymax = upper)

(p <- ggplot(plot.dta, aes(rep.type, effect)) +
    facet_wrap(~city.type) +
    geom_errorbar(limits, width = 0.1) +
    geom_line() +
    theme_bw(base_family = "serif", base_size = 12) +
    theme(panel.spacing = unit(2, "lines"),
          plot.margin = unit(c(1,2,1,1), "cm"),
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white")) +
    labs(y = "Probability of Lobbying", x = "Average Delegation Conservatism Score"))

